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Abstract 

We present investigations of the combined effects of Debye-Hiickel reput- 



ed 
H ■ 

sive and overlaping Debye spheres attractive interaction potentials around 
charged dust particles on collective modes, phase separation and ordered 
structures in a strongly coupled dusty plasma. We obtain static and dy- 
namical information via Molecular Dynamics simulations in the liquid and 
crystallized phases and identify the onset of an instability in the transverse 
mode, by using lattice summation method. The results are useful for under- 
standing the origin of coagulation/agglomeration of charged dust particles 
and the formation of ordered dust structures in low-temperature laboratory 
and space plasmas. 
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1. Introduction 



Dusty plasmas are composed of electrons, ions and dust particles of dif- 
ferent sizes (ranging from nanometer to micrometer). Dust particles, which 
are negatively charged by absorbing electrons from the background plasma, 
respond to electrostatic and electromagnetic fields and become active com- 
ponents of the dusty plasma system, especially when their dynamics comes 
into the picture at kinetic levels l|, |2( . Dusty plasmas are ubiquitous in many 
physical environments, viz. in interstellar media and nebulas, in our solar 
sysem (e.g. Saturn's ring system, cometary tails and comas, in in the Earth's 
mesosphere and ionosphere), in industrial and laboratory processing plasmas 
for high-technological applications (also in medicine), as well in laboratory 
gas discharges and magnetically confined fusion reactors. In low-temperature 
laboratory dusty plasmas on ground and on International Space Station un- 
der microgravity, dust is intentionally introduced in plasmas for fundamental 
studies of collective phenomena at kinetic level (e.g. the discovey of the dust 



acoustic wave {3]), involving an ensemble of charged dust grains. In astro- 
physics, the coagulation of charged dust particulates plays a significant role 
in the formation of large scale structures (e.g. plenetesimals) . There are 
forthcoming space missions for exploring the composition and dynamics of 
dust particles in the vicinity of Sun, as well on Mars and Moon. 

To describe the interaction between the charged dust particles in dusty 
plasma systems, the Debye-Hiickel (DH) potential <j) DH (r) has been widely 
adopted as a model potential around a chargd dust grain that has a well 
defined surface potential. The DH potential expresses that the Coulomb 



repulsion between the charged gains of similar polarity is screened by the 
plasma environment, where electrons and ions represent a polarizable back- 
ground. Consequently, the 1/r Coulomb potential is exponentially shielded 
over a characteristics distance, the dusty plasma Debye radius |l| Xn, and 
we have 

0^(r) = -^ eXp( - r/Ap) , (1) 
Arce r 

where Q is the charge of the dust particles and e is the permittivity of free 
space. 

At short inter-particle distances, the DH interaction potential has to be 
modified, since the Debye spheres of the neighboring grains would overlap 
{4]. From classical electrostatics, one can then derive an additional term 
of the interparticle potential that expresses the interaction between these 
overlapping Debye spheres (ODS) |4j,|5j. The resulting (total) potential in 
this case becomes 



r DS (r)- Q 



1 1 



Wr 2A D J «P(-'A>>. P) 
which exhibits an attractive force beteen charged dust grains of the same 
polarity at distances beyond r = v3Ai> 

In this Letter, we investigate the effect of this peculiarity of the attrac- 
tive potential on the structure of two-dimensional dust-plasma systems in 
the crystallized phase and in the strongly coupled liquid phase, where the 
potential energy dominates over the kinetic energy, i.e. r ^> 1, where V is 
the coupling parameter, viz. 
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r = ^ 2 (3) 

Ane Q ak B T' y ' 

Here kg is the Boltzmann constant, T the dust temperature, and a = 

{\/im) l l 2 the Wigner-Seitz radius (that characterizes the inter-dust particle 
spacing), with n being the dust particle number density. We also introduce 
the dimensionless screening parameter 

k = f . (4) 

D n 

In a previous theoreical investigation 6|, the attractive term in (fi ODS (r) 
has been shown to cause instability of transverse dust lattice oscillations. 
Considering the nearest neighbor interactions only, the onset of instability 
was identified at b/Xo > v3+ 1, where b is the lattice constant. For the 
triangular lattice b rs 1.905a, the critical screening value in terms of the 
inter-grain spacing becomes k c = a/Xo ~ 1.43. 

As regards to the strongly coupled liquid phase, Molecular Dynamics 
simulations have confirmed that <ft ODS (r) gives rise to an uneven spatial dis- 
tribution (agglomeration) of dust particles [7|. 

2. Methods 

In our numerical work, we use a standard Molecular Dynamics simulation 
method with Langevin dynamics js, 9]. The system is two-dimensional and 
we use periodic boundary conditions. We carry out computations with both 
potentials, <ft DH (r) and <p ODS {r). The simulations are initialized either by a 
spatially random particle configuration (in the studies of the liquid phase, 
with iV = 10,000 dust particles) or by a lattice configuration (in the studies 
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of the crystallized phase, with N = 9180 dust particles). The side lengths of 
the simulation box are chosen to accommodate a perfect triangular lattice. 
Initial dust particle velocities are sampled from a Maxwellian distribution 
of the temperature T , corresponding to the specified value of the coupling 
parameter T. The interaction between the dust particles is limited to dis- 
tances smaller than a cutoff radius, r*, which is chosen to be large enough so 
that interparticle forces at r > r* become very small due to the exponential 
term in the interaction potential. The neighbors of the dust particles within 
r < r* are searched using the chaining mesh technique. 

The equation of motion includes terms representing the frictional damping 
and random forces, which express random collisions of gas molecules with the 
charged dust particles. We have 

dv- N 

m —T~ = ~Q V0ij — vmw + R, (5) 

where <pij is the interaction potential between dust particles i and j, v is 
the frictional damping, and R is a random force, see e.g. {§], 9]. In all 
our simulations, we use a small friction value, G = v/ojq = 0.01, where 
ujq = a/ nQ 2 /2e ma is the nominal 2D dust plasma frequency, with m being 
the dust particle mass. We introduce dimensionless distance as f = r/a and 
the wave number as k = ka. 

The pair correlation function and the static structure function are derived 
in the simulations in the standard manner. Information about the (thermally 
excited) collective modes and their dispersion is obtained from the Fourier 
analysis of the correlation spectra of the density fluctuations, 
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yielding the dynamical structure function as [10J: 



(6) 



S(k, (J) 



(7) 



where At is the length of data recording period and p(k,u>) = J-\p(k,t)] is 
the Fourier transform of (jSJ). 

Similarly, the spectra of the longitudinal and transverse current fluctu- 
ations, L(k,u) and T(k,u) can be obtained from Fourier analysis of the 
microscopic quantities, respectively, 



where Xj and Vj are the position and velocity of the j-th particle. Here 
we assume that k is directed along the x axis (the system is isotropic) and 
accordingly omit the vector notation of the wave number. This calculation 
allows the derivation of the spectra for a series of wave numbers, which 
are multiples of k miritX ^ = 2n/L x (y), where L x ^ is the edge length of the 
simulation box in the x (or y) direction. The collective modes are identified 
as peaks in the fluctuation spectra. 

The phonon dispersion in the lattice configuration is calculated in terms 
of the lattice dynamical matrix defined as 





(8) 
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c,M 



Q 



(9) 



m 




with a summation over all the lattice sites j, keeping i fixed (r^ = 0). The 
lattice dynamical matrix reflects the symmetry of the underlying lattice. The 
diagonalization of is possible in the coordinate system of the eigenvectors, 
whose orientations, in general, do not coincide either with the direction of 
k or with the crystallographic axes. To find the eigenmodes, we follow the 
traditional method of solving the secular equation 



3. Results 

In figured] we present an overview of the main static characteristics of our 
dusty plasma system: the four panels display the pair correlation function 
[g(r)} and the static structure function for the 0.5 < k < 3 domain 

of the screening parameter, for both types of potentials. The coupling pa- 
rameter is fixed at T= 120; at this value the system is in the liquid phase 
for the whole range of k values covered. Panels (a) and (b) correspond to 
4> DH , while (c) and (d) correspond to cfP DS . At small screening values, viz. 
k < 1, the g(r) and S(k) functions are similar for both potentials, showing 
slightly less organized structures (indicated by smaller amplitude the peaks) 
in the case of (fP DS . In the case of (f> , we observe major changes in both 
static quantities in the vicinity of the screening value k w 1.3. In g(r) a 



u; 2 -CV(k)|| =0. 



(10) 
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Figure 1: Static characteristics of the 2D dusty plasma obtained with <f) DH (top row) and 
with cjP DS (bottom row). (a,c) pair correlation function, g(r/a), (b,d) static structure 
function S(k). T = 120. 

strong peak develops at small separations, indicating particle agglomeration 
and S(k) suddenly increases at small wave numbers, indicating development 
of large-scale structures. 

Figure |2] shows a comparison between systems characterized by <ft DH and 
4> , at T = 100 and k — 2. At this screening strength, the system with 
(j) already exhibits agglomeration, as the particle snapshot and the corre- 
sponding pair correlation function (in panels (e) and (f) of figure [2]) illustrate. 
A notable difference in the L(k, u) spectrum is the significant increase of the 
maximum mode frequency from u/uq ~ 0.5 (at <p DH ) to wl.l (at <f) ODS ), 
and an increase of the wave number where this frequency is found, from k m 
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Figure 2: MD simulation results obtained with the potentials cj) DH (top row) and <fr ODS 
(bottom row). (a,e) Snapshots of particle positions, (b,f) pair correlation functions, (c,g) 
spectra of longitudinal current fluctuations, and (d,h) spectra of transverse current fluc- 
tuations. T = 100, k = 2. 

1.9 to ~ 3.1. As the transverse mode is quite weak in the liquid phase, it is 
more difficult to characterize quantitatively the frequency of this mode, but 
an increase of the frequency is observable, too. 

Next, we compare the results of our MD simulations with those obtained 
from lattice summation. To be able to do that, we have initialized the simula- 
tions with a lattice configuration and set a very high coupling value, T = 10 4 . 
Figure |3] shows the results for three cases: (i) <p DH with k = 1, (ii) <f) ODS with 
k = 1.47, and (iii) with k = 1.48 The lattice is stable in the first two of 
these three cases, while the third case belongs to the unstable domain with 
. The MD and lattice calculation results are in very good agreement 
for the stable cases. The figure exhibits that the L and T modes are super- 
imposed, at propagation directions of a = 0° (the direction of the nearest 
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Figure 3: Sum of the current fluctuation spectra, L(k, u) + T(k, w), derived from the MD 
simulations (color maps) and dispersion relations of the collective excitations obtained 
from zero-temperature lattice summation (solid lines: black: longitudinal modes, white: 
transverse modes). Upper row (a,c,e): propagation direction: a — 0°, bottom row (b,d,f): 
a = 30°. (a,b) 4> DH with k = 1, (c,d) (f> ODS with k = 1.47, (e,f) (j> ODS with k = 1.48. All 
MD runs have been initialized with lattice configuration of particles, at T — 10 4 . For the 
color maps the scale is logarithmic and covers three orders of magnitude. 
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Figure 4: Results of lattice summations for the eigenfrequencies of the two modes, for 
<p . Red color denotes quasi- longitudinal and blue color denotes quasi-transverse modes; 
note that the polarization of the modes is purely longitudinal or transverse in the principal 
directions (a = 0° and a = 30°) only. The panels show the frequencies in propagation 
directions between a = 0° and a = 30° with 1° resolution, (a) Stable domain, k = 1.49. 
(b) Zoomed part of (a) showing the square of the mode frequencies at small wave numbers. 
In this domain, all dispersion curves for a = 0° . . . 30° overlap, (c) Same as (b), but for 
k = 1.50; note the splitting of the dispersion curves and the appearance of negative (w/wo) 2 
values, indicative of an instability. The arrow shows increasing values of a. 

neighbor) and a = 30°. In these propagation directions, the polarization of 
the modes becomes purely longitudinal and transverse, in other directions 
the polarizations mix. At k = 1.49 complicated structures appear in the 
fluctuation spectra that are difficult to interpret. 

Figure H] displays the results of lattice summations for the eigenfrequencies 
of the two modes, for <p ODS , characterizing the interaction of dust particles 
with overlapping Debye spheres. Figure IU(a) shows the dispersion relation 
of the modes with different propagation directions, for k = 1.49, which be- 
longs to the stable regime of the lattice at T = 0. The two branches of 
the curves correspond to "quasi-longitudinal" and "quasi-transverse" modes. 
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The polarization of the modes is usually mixed, except for the principal di- 
rections (a = 0° and a = 30°), where the polarization is purely longitudinal 
or transverse. Figure IU(b) exhibits (u/ojq) 2 values for small wave numbers 
at k — 1.49. In this domain, all dispersion curves for a = 0° . . . 30° overlap. 
Figure 0(c) shows also a zoomed part of the mode dispersion relations, but 
for k = 1.50; at this screening value negative (u/uo) 2 values appear that 
correspond to an instability of the lattice. 

The critical value of the screening parameter was found to be k c = 1.4997 
in our lattice calculations (that add contributions to the dynamical matrix 
up to the 126th neighbors of the particles). The theoretical calculations 
taking into account the nearest neighbor interactions only resulted in k c « 
1.43. In our MD simulations carried out for T = 10 4 we found that the 
lattice disintegrates between k = 1.47 and k = 1.48, while in the fluid phase, 
as shown in Figured] for T = 120, an instability occurs at k c rj 1.3. These 
changes indicate the effect of the temperature on the development of the 
instability. 

4. Summary 

Summing up, we have investigated the effect of an attractive term of 
the interaction potential in addition to the DH potential, conventionally 
used to model dusty plasmas, on the static and dynamical characteristics 
of many-body dust-plasma systems. Molecular dynamics simulations (with 
a Langevin approach) carried out in the liquid phase showed that the onset 
of particle agglomeration at screening values exceeding a critical value. The 
particle agglomeration was accompanied with a significant increase of the 
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frequencies of the longitudinal and transverse collective excitations. In the 
crystallized phase, the MD simulations revealed that the hexagonal lattice 
becomes unstable at high screening. These observations have been confirmed 
by harmonic lattice summations in the zero temperature limit, which showed 
the development of an instability in the (quasi-)transverse modes at the crit- 
ical screening value of k c m 1.4997. In conclusion, we stress that the results 
of the present investigation are useful for understanding the origin of agglo- 
maration/coagulation of charged dust particles and the formation of ordered 
dust structures in low-temparature laboratory and space dusty plasmas 2|. 

This work was supported by the Hungarian Fund for Scientific Research 
through grants K77653, K105476, NN103150, and by the HELIOS project 
(ELI 09-1-2010-0010). 
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